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ABSTRACT 

To date, most simulations of the chemistry in protoplanetary discs have used 1+1D 
or 2D axisymmetric G-disc models to determine chemical compositions within young 
systems. This assumption is inappropriate for non-axisymmetric, gravitationally un¬ 
stable discs, which may be a significant stage in early protoplanetary disc evolution. 
Using 3D radiative hydrodynamics, we have modelled the physical and chemical evo¬ 
lution of a 0.17 M 0 self-gravitating disc over a period of 2000 yr. The 0.8 M 0 central 
protostar is likely to evolve into a solar-like star, and hence this Class 0 or early Class 
I young stellar object may be analogous to our early Solar System. Shocks driven by 
gravitational instabilities enhance the desorption rates, which dominate the changes 
in gas-phase fractional abundances for most species. We find that at the end of the 
simulation, a number of species distinctly trace the spiral structure of our relatively 
low-mass disc, particularly CN. We compare our simulation to that of a more mas¬ 
sive disc, and conclude that mass differences between gravitationally unstable discs 
may not have a strong impact on the chemical composition. We find that over the 
duration of our simulation, successive shock heating has a permanent effect on the 
abundances of HNO, CN and NH3, which may have significant implications for both 
simulations and observations. We also find that HCO + may be a useful tracer of disc 
mass. We conclude that gravitational instabilities induced in lower mass discs can sig¬ 
nificantly, and permanently, affect the chemical evolution, and that observations with 
high-resolution instruments such as ALMA offer a promising means of characterising 
gravitational instabilities in protosolar discs. 

Key words: stars: pre-main-sequence, stars: circumstellar matter, protoplanetary 
discs, astrochemistry 


1 INTRODUCTION 

In the very early stages of stellar formation, the surrounding 
protoplanetary disc can contain a mass comparable to that 
of the central protostar (e.g. Machida & Matsumoto 2011). 
As a result, self-gravity within the disc is significant and the 
disc can become unstable if gravitational forces overcome 
pressure and shear forces; a criterion characterised by the 
Toomre parameter (Toomre 1964). Gravitational instabili¬ 
ties (GIs) formed in the disc via this process produce spi- 
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ral density waves that efficiently drive angular momentum 
transport outwards and mass accretion inwards. As these 
spiral waves grow, they can produce shocks that heat the 
disc material locally (e.g. Harker & Desch 2002; Boley & 
Durisen 2008; Bae et al. 2014). This episodic heating can en¬ 
hance the rates of chemical reactions, such as desorption of 
volatiles and reactions with high activation energies, which 
may significantly affect the chemical composition of discs. 

Many authors have investigated the chemical proper¬ 
ties and evolution of protoplanetary discs (see Henning & 
Semenov 2013, Table 3 for an extensive review). However, 
most simulations have featured a passive or weakly vis- 
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cous a-disc with 2D axisymmetric structure and local mass 
transport. Hence, this prescription is likely inappropriate 
for early phases of disc evolution if GIs are present driving 
non-axisymmetric structure. To date, only one study has 
explored the effect of GIs and subsequent shock-heating on 
disc chemistry. 

Ilee et al. (2011), hereafter 12011, simulated a massive 
protoplanetary disc throughout an evolutionary phase that 
may be applicable to an FU Orionis outburst. The system 
represents a Class 0/early Class I object and the relatively 
massive (0.39 M 0 ) disc is appropriate for a protostar that 
will likely become an F star (« 1.4 M©) at the end of the ac¬ 
cretion phase. The authors found that shock-heating within 
the disc induced temperature peaks that desorbed various 
chemical species from the surface of dust grains and en¬ 
hanced the rates of reactions that were otherwise not en¬ 
ergetically favourable, therefore increasing some gas-phase 
abundances. In contrast to simulations of more evolved, 
less dynamic systems, they found that the midplane is the 
hottest region and hence rich in chemistry. All species that 
they consider are tracers of the spiral shocks, and Douglas 
et al. (2013) showed that, encouragingly, these gas-phase 
tracers of disc dynamics should be detectable using ALMA 
at a distance of 100 pc. Furthermore, Dipierro et al. (2014) 
have recently shown that ALMA can readily detect spiral 
structure in continuum emission, within Class 0/1 systems 
located at distances comparable to TW Hydrae, Taurus- 
Auriga and Ophiucus star-forming regions. Hence this will 
allow the empirical assessment of whether young protoplan¬ 
etary discs around intermediate mass stars are indeed grav¬ 
itationally unstable. 

Simulations of gravitationally unstable discs can be 
used to determine which observable features are the most 
promising for detection, should the instabilities be present. 
Observations can then test and constrain models for the im¬ 
portant dynamical and chemical evolution of young proto¬ 
planetary discs. If gravitational instabilities are not detected 
in real objects, then the model parameters can be refined un¬ 
til GIs can be satisfactorily excluded as a possibility in young 
systems, or the simulations can be improved to identify ac¬ 
tual observable signatures. The motivation for this paper is 
to provide a first step in this direction for a protosolar disc, 
i.e. a young disc surrounding a protostar that will become a 
Solar-like star at the end of the accretion phase. 

In Section 2 we outline the physical and chemical mod¬ 
els we have used and how they differ from those used in 
the simulation of the more massive disc. Section 3 contains 
results for time-dependent fractional gas-phase abundances 
determined from individual fluid elements and column den¬ 
sity maps of the entire disc in order to assess the differences 
in composition between the more massive disc, the lower 
mass disc at different times, and discs featured in other pub¬ 
lished works. Finally, in Section 4 we present conclusions 
based on our results and discuss future research. 


2 DISC MODELS 

2.1 Dynamical simulation 

A radiation hydrodynamics simulation of a 0.17M© proto¬ 
planetary disc is run over a period of approximately 2000 yr 


and used as the input into our chemical model. The disc 
surrounds a central protostar of 0.8 M 0 , which we envisage 
will evolve into a Solar-like star. The majority of the mass 
in the disc initially ranges between r « 6-41 au from the 
central protostar, but spreads partially during its evolution 
as a consequence of angular momentum transport and ac¬ 
cretion, resulting in a disc that spans r « 5-54 au at the end 
of evolution. The system represents a Class 0 or early Class 
I object that can potentially be compared to the early Solar 
System. 

The physical simulation is performed using 
CHYMERA (Boley 2007), which solves the equations 
of hydrodynamics on a regularly spaced cylindrical grid. 
Outflow boundaries are used at the inner, outer, and top 
grid edges and mirror symmetry is assumed along the 
midplane of the disc. Self-gravity is calculated directly 
through cyclic reduction, which is a numerical method for 
solving large linear systems by repeatedly splitting the 
problem, with the boundary potential determined by a 
spherical harmonics expansion. An indirect potential is 
used for capturing the star-disc interactions (e.g. Michael 
& Durisen 2010) and the Boley et al. (2007a) equation of 
state is used for a fixed H 2 ortho-para ratio of 3:1 in a 
Solar mixture gas. We use the Boley et al. (2007b) radiative 
transfer algorithm, which combines flux-limited diffusion 
with ray tracing in the vertical direction. The radiative 
routine calculates a separate radiative time step, and 
subcycles in the event that the hydrodynamics time step 
becomes much larger than the radiative step. A maximum 
iteration of 8 sub cycles is used. If this limit is reached, then 
one last step is used to synchronize the hydrodynamics and 
radiative times. To ensure stability, energy is not allowed to 
change in any one cell by more than 10 per cent (see Boley 
& Durisen 2006); such limiters are typically only necessary 
in very low-density regions. The disc is assumed to be 
irradiated by the protostar, with a background temperature 
profile Ti rr 4 = (150 K(r/au) 0 ' 5 ) 4 + (3K) 4 , where r is radius 
from the central protostar. 

The disc is initialized with a flat Q ~ 1.2 profile 
(Toomre 1964) and a temperature profile following T oc 
r -0 ' 5 , which gives a surface density profile E oc r _1 ‘ 75 for 
a Keplerian disc. The disc is seeded with a random density 
perturbation to promote the growth of gravitational insta¬ 
bilities. The initial Q value is low enough to guarantee a 
very rapid onset of spiral structure, but not so low as to 
overshoot the instability regime severely. While the 12011 
calculations were performed in the context of a massive disc 
undergoing an outburst of GIs, the simulations here are in¬ 
tended to explore a more protracted phase of instability in 
a lower mass disc. As such, we evolve the new disc past 
the outburst phase of unstable discs (Mejfa et al. 2005) be¬ 
fore we begin tracing fluid elements. After about 1290 yr 
of evolution (about 5 orbits at 40 au), 2000 fluid elements 
are placed randomly in the disc, but weighted by mass to 
reflect the actual distribution of material. This is twice as 
many fluid elements as used in 12011. The fluid elements 
are evolved by interpolating the gas flow from the surround¬ 
ing cells to a given tracer’s position and integrated directly. 
Local disc conditions such as pressure, density, and temper¬ 
ature are also interpolated to the fluid element’s position, 
allowing us to trace the thermodynamic history of the gas. 
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Figure 1 . Nuclei column density at t = t 0rp for the 0.17Mq 
(left) and 0.39 Mq (right) discs as viewed from above. The more 
massive disc (right) undergoes a violent instability, while the lower 
mass disc (left) shows a more protracted phase of instability. 



Figure 2. As in Figure 1, but showing the maximum line-of-sight 
temperature. 


Of the 2000 fluid elements evolved in the simulation, nine 
are lost through a grid boundary and are not included in 
the following analysis. 


2.2 Comparison with the more massive disc 

The more massive disc featured in 12011 was simulated for 
389 yr, and because of this short duration the hydrodynam- 
ical results were repeated for ten cycles, i.e. using the final 
abundances after 389 yr as initial abundances for the next 
run, in order to achieve a well mixed chemical distribution. 
This introduced unrealistic discontinuities into their chem¬ 
ical calculations. Here, we run the new simulation continu¬ 
ously for 2047 yr with fluid element tracers present, which 
means we can follow the behaviour of the lower mass disc 
over a longer period without artificial cycling, allowing for 
a more self-consistent analysis. 

The more massive disc was initialised with a Q ~ 1 pro¬ 
file for most of the disc and a TM — 140(r/au)°' 5 + 10 K ra¬ 
diation field, requiring a surface density profile E oc r -1 ' 75 . 
Although these parameter initialisations are similar to the 
lower mass disc setup, these discs are not directly compa¬ 
rable because they trace different potential phases of proto¬ 
planetary disc evolution; the higher mass disc is more appli¬ 
cable to a violent burst in activity, whereas the lower mass 
disc undergoes a more quiescent protracted period of evolu¬ 


tion. However, we can compare the discs after a consistent 
number of outer rotational periods (ORPs), which defines 
the time since the inclusion of tracers at a particular radii. 
At 30 au the ORP is approximately 140 yr in the more mas¬ 
sive disc and 165 yr in the less massive disc. We find that 
the more massive disc has completed 2.7 ORPs at 30 au by 
the end of the simulation, t = 389 yr. The lower mass disc 
completes the same number of ORPs at the same radius af¬ 
ter approximately t = 440 yr. Hence, the discs are roughly 
comparable dynamically at these respective times, which we 
denote as t = t orp for both discs, and we use this as the basis 
for our comparisons. 

Figures 1 and 2 show the column density of total nuclei 
and maximum line-of-sight temperature, respectively, within 
the featured discs at t — t 0rp , as viewed from above. At 
this time the maximum temperature and number density 
within the more massive disc are 400 K and 8 x 10 13 cm -3 , 
respectively, which is 1.2 times hotter and 1.4 times denser 
than within the less massive disc. Both figures show that 
prominent spiral features have formed in each disc but the 
more massive disc is significantly more flocculent due to the 
stronger gravitational instabilities and evolutionary phase. 

Figure 3 shows slices of the disc interior along the x = 
Oau and y = Oau axes within the less massive and more 
massive discs at t — t 0rp . Overall the interior structures 
are similar; both discs show a flared morphology, with the 
hottest and densest material located in a ring-like structure 
approximately 10 au from the protostar, and possess uneven 
surfaces resulting from the shocks (Boley & Durisen 2006). 
However, the more massive disc is also clearly more extended 
and features enhanced temperatures at large radii, in con¬ 
trast to the lower mass disc. Moreover, there is a higher 
ambient temperature because of the more massive central 
protostar and disc. These differences in outer disc temper¬ 
ature structure in particular may have a prominent effect 
on the chemistry as more energy is available for chemical 
processes and reactions at larger radii. 

The number densities and temperatures of fluid ele¬ 
ments tracked throughout the disc simulation are used as 
the input for the chemistry model. Figure 4 shows the tra¬ 
jectories of a fluid element in the outer regions of each disc, 
which is only intended to be illustrative of the types of tra¬ 
jectories that fluid elements can follow. Both fluid elements 
are initially roughly 30 au from the centre and orbit on a 
trajectory that first radially expands and then contracts. 
We identify the collisional shocks experienced by each fluid 
element from the concurrent peaks of temperature and num¬ 
ber density, seen in Figure 5, and pressure (not shown). The 
shocks cause the fluid elements to be displaced vertically 
upwards in the discs before falling back towards the mid¬ 
plane. This is an effect seen in gravitationally unstable disc 
simulations as the strong spiral waves in a vertically strat¬ 
ified disc are both shocks and hydraulic jumps (e.g. Boley 
& Durisen 2006). The velocities of the shock waves in the 
disc are relatively low (a usual Mach number is 2, implying 
a shock velocity of a few km s -1 ). Such low Mach numbers 
are due, in part, to the shallow pitch angles of the spiral 
arms, which significantly reduce the speed of a fluid element 
normal to the spiral shock (see Boley & Durisen 2008, for a 
detailed discussion). 
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Figure 3. Nuclei number densities and temperatures within slices of the disc interior at t = t 0 rp for the 0.17 Mq (top four panels) and 
0.39 Mq (bottom four panels) discs. Slices are taken at the specified location and d denotes the distance from the disc centre along the 
appropriate orthogonal axis. 




x [au] 




Figure 4. Positions and corresponding temperatures of fluid elements extracted from the 0.17 M© (top) and 0.39 M© (bottom) discs 
that demonstrate similar orbital behaviour within 2.7 0RPs at 30 au in each disc. The open circles denote t = 0 and the closed circles 
denote t = t Q rp- 


2.3 Evolution of the lower mass disc 

Figures 6 and 7 show the nuclei column density and max¬ 
imum line-of-sight temperature, respectively, within the 
0.17 M© disc at t — t Qrp (t — 440 yr) and at the end of 
the simulation t — (t = 2043yr), as viewed from above. 
Over the duration of the simulation the average mass flux 
through the disc is 10 -6 M©yr -1 . By the end of the simula¬ 


tion the disc is clearly more flocculent as the gravitational 
instabilities have driven progressively denser spiral waves. 
A video of the lower mass disc evolution demonstrates this 
and is provided as online supplementary material. 

Figure 8 shows the trajectories of two fluid elements lo¬ 
cated within the inner and outer regions of the lower mass 
disc, throughout the entire duration of the simulation. This 
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Figure 5. Temperature and number density histories of fluid ele¬ 
ments extracted from the 0.17M© and 0.39 M© discs that follow 
similar trajectories within 2.7 ORPs at 30 au in each disc (see 
Figure 4). 


3 

03 


log N(n) [cm- 2 ] 



Figure 6. Nuclei column density at t = t 0rp (left) and t = 
(right) for the 0.17 M© disc as viewed from above. As the disc 
evolves dynamically, the spiral structure becomes much more floc- 
culent. 
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Figure 7. As in Figure 6, but showing maximum line-of-sight 
temperature. The maximum temperature reached is approxi¬ 
mately 600 K. 


behaviour is also shown in the movie provided as online 
supplementary material (blue trails). The temperature and 
number density histories of these fluid elements are shown 
in Figure 9, and the outer disc temperature and number 
density evolution displayed is for the same fluid element as 
featured in Figures 4 and 5. In the outer disc, the temper¬ 
ature profile is determined primarily by the stellar irradia¬ 
tion, with excursions due to spiral shock heating. As a result, 
over the duration of the entire simulation, the temperature 
and number density of the outer disc fluid element appear 
periodic. This period is approximately 100 yr, which means 
that the spiral waves propagate faster than the outer disc 
rotates; the average corotation of the spiral structure ap¬ 
pears to be at about 17 au based on the average mass flux 
profile. The inner disc fluid element, however, exhibits sig¬ 
nificantly different behaviour. The much faster rotation of 
this fluid element, roughly 10 yr per orbit, results in much 
more frequent shock encounters, which produces the signif¬ 
icantly faster variations in temperature and number den¬ 
sity. Moreover, there is a substantial increase in the average 
temperature and number density for the inner disc fluid el¬ 
ement after about t = 400 yr. This is because the gravita¬ 
tional instabilities drive mass transport and shock heating. 
Furthermore, within « 20 au, the radiative cooling becomes 
very inefficient as the densest, inner disc regions are more 
optically thick than the outer disc. In the inner disc fluid 
element, the rapid variations in temperature and number 
density, approximately every 10 yr, reflect passage through 
spiral structure. The large-scale changes, however, spanning 
over 100 yr or more, reflect different regions of the inner disc 
that are in rough pressure equilibrium, i.e. cool and dense 
regions, or hot and rarefied regions. 


2.4 Chemical Model 

We use the chemical evolution code developed in 12011, 
which consists of a network of 125 species and 1334 re¬ 
actions. These reactions are selected from a subset of 
UMIST95 (Millar et al. 1997), whose rate coefficients have 
been updated using KIDA 1 , and consist of many processes: 
positive ion-neutral reactions, ionisation by cosmic rays, 
charge transfer, proton transfer, hydrogen abstraction, ra¬ 
diative and dissociative recombination with electrons, ra¬ 
diative association, neutral exchange, photodissociation and 
photoionisation, recombination of ions with negative grains, 
and adsorption and desorption. The code interpolates tem¬ 
perature and number densities extracted from the fluid el¬ 
ements to a higher resolution chemical time-scale using cu¬ 
bic spline fits. The DVODE package (Brown et al. 1989) is 
then used to integrate the rate equations for each considered 
species, outputting fractional gas-phase abundances, Xi — 
rii/n, where n — tih + riHe + riz (see 12011, Equation 1). 

While the chemical model includes various photochemi¬ 
cal reactions, we perform the calculations under the assump¬ 
tion that the disc environment of Class 0 or early Class I 
objects is well shielded from sources of stellar and inter¬ 
stellar radiation, and therefore take Ay = 100 mag. The 
combination of envelope infall onto the disc and powerful 


1 http://kida.obs.u-bordeauxl.fr/ 
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Figure 8. Orbital positions and corresponding temperatures of fluid elements extracted from the inner (top) and outer (bottom) regions 
of the 0.17 M© disc. The open circles denote t = 0 and the closed circles denote t = tfi n . 


outflows (see, e.g., Machida V Hosokawa 2013) may dras¬ 
tically limit the illumination of the upper layers of young 
discs by the protostellar and interstellar UV field. However, 
we allow photochemical reactions induced by cosmic rays, 
the ionisation rate of which we take to be £ = 10 -17 s -1 . 

The chemical model assumes that simple surface chem¬ 
istry (i.e. hydrogenation) occurs immediately as species are 
adsorbed on to the grain surfaces. In order to simulate any 
further potential effects of grain surface reactions, our ini¬ 
tial species abundances are taken to be representative of ob¬ 
servations of cometary ice abundances (see Ehrenfreund & 
Charnley 2000, Table 2) and are given in Table 1. Although 
there is some consistency between cometary and interstel¬ 
lar ices, it is not known whether comets undergo significant 
chemical processing, and thus, whether they accurately pre¬ 
serve the initial elemental compositions of young systems is 
not clear (see Caselli V Ceccarelli 2012, Chapter 7). There¬ 
fore, in the future we shall include a detailed treatment of 
grain surface processes in our chemical model in order to 
quantify their effect. 

As mentioned previously, the velocities of the shock 
waves in the disc are relatively low (a few kms -1 ). Thus, 
we do not expect these shocks to cause effects such as the 
disruption of dust grain cores or sputtering, which typically 
require shock speeds of at least 25kms -1 (e.g. Caselli et al. 
1997; Van Loo et al. 2013). Rather, we focus on the effects 
of the shocks on the thermal desorption of the ices from the 
icy mantles of dust grains. 


2.4.1 Modifications 

The chemical evolution code that we use was first presented 
in the 12011 paper. However, we make some alterations in or¬ 
der to improve its efficiency and success rate. Originally the 
code used a logarithmic time step in order to capture the 
fastest initial chemical reactions at a sufficient resolution. 



Figure 9. Temperature and number density of fluid elements 
extracted from the outer and inner regions of the 0.17M© disc 
(see Figure 8). The rapid variations in temperature and number 
density of are a consequence of collisional shocks, which are much 
less frequent in the outer disc. The large-scale changes seen only 
in the inner disc are a result of the fluid element passing through 
cooler and denser, or more rarefied and hotter regions of the inner 
disc that are in rough pressure equilibrium. 

This was necessary because the fluid element temperature 
and density histories in the more massive disc were cycled to 
compensate for the short time-scale of the physical disc evo¬ 
lution. This introduced artificial numerical discontinuities 
into the chemical analysis that caused the chemical integra¬ 
tor to fail for a selection of fluid elements. The less massive 
disc is simulated for a considerably longer time, so we have 
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Table 1. Initial fractional abundances X(i) = rii/n. Note that 
a(b) represents a x 10 b . 


Species 

Abundance 

He 

1.00(-1) 

H 2 CO 

1.83 (- 6 ) 

co 2 

3 . 67 (- 5 ) 

CO 

3 . 66 (- 5 ) 

HCN 

4.59 (- 7 ) 

ch 4 

1 . 10 (- 6 ) 

HNC 

7 . 34 (- 8 ) 

S 

1.62(-5) 

nh 3 

3 . 30 (- 6 ) 

h 2 s 

2 . 75 (- 6 ) 

so 

1 . 47 (- 6 ) 

h 2 o 

1.83 (- 4 ) 

so 2 

1 . 84 (- 7 ) 

OCS 

3 . 30 (- 6 ) 


no need to cycle the chemical abundances. Therefore, we 
changed the time step to be additive, with a time step of 
1000 seconds so that fast reactions were followed with suffi¬ 
cient resolution throughout the entire disc evolution. The 
resulting abundances were recorded every 0.05 yr. Recur¬ 
sion was also incorporated into the chemical code so that 
a failed integration step was reattempted in a dichotomy 
paradox fashion, i.e. continually halving the time-step, until 
the integration succeeded or the step size became zero and 
the run failed. These two adaptations made the code much 
more robust against rapid changes in temperature and den¬ 
sity, which a significant number of inner disc fluid elements 
possess. 


3 CHEMICAL RESULTS 

Understanding how gravitational instabilities affect the 
chemical evolution of the lower mass disc, and how this dif¬ 
fers from the more massive disc, is necessary for understand¬ 
ing the potential diversity and observational signatures of 
such discs. This allows us to determine what tests are the 
most likely to detect or exclude active instability in any 
given Class 0 or I source. 

Figure 10 shows the fractional gas-phase abundance for 
20 species (CO, C0 2 , SO, S0 2 , NH 3 , H 2 0, H 2 S, H 2 CO, 
OCS, 0 2> HCO, HCO+, HNO, CS, HCS, HCS+, HCN, HNC, 
CN and OCN) as a function of t 0rp for both discs. The abun- 
dances are shown for the same fluid elements featured in 
Figures 4 and 5, which as discussed orbit at approximately 
30au. We have included three species, C0 2 , CN and H 2 CO, 
in addition to those that are featured in the 12011 paper 
because these are significant species that have recently been 
detected in protoplanetary discs (e.g. Chapillon et al. 2012 ; 
Qi et al. 2013; Bast et al. 2013). 

The gas-phase fractional abundances of the species in 
the upper right panels (CO, SO, OCS, NH 3 , H 2 0 and H 2 S) 
are only significantly affected by adsorption and desorption 
processes because their total abundance, comprised of gas- 
phase, X(i), and grain abundance, N g (i)rj , is approximately 
constant during the disc evolution. Here, 77 is the ratio of 
the number density of dust grains to the number density of 


Table 2. Maximum fractional gas-phase abundances within the 
0.39 Mq disc at t — t 0 rp, 0.17 M© disc at t = t 0rp and O.17M0 
disc at t = . 


Species Maximum log X(i) 


i 

0.39 M© 

t = torp 

0.17M© 

t = torp 

0.17 M© 
t = ffin 

CO 

-4.4 

-4.4 

-4.4 

co 2 

-4.4 

-4.4 

-4.4 

so 

-5.8 

-5.8 

-5.8 

so 2 

-6.5 

- 6.6 

-6.3 

nh 3 

-5.5 

-5.5 

-5.5 

h 2 o 

-3.7 

-3.7 

-3.7 

h 2 co 

-5.7 

-5.7 

-5.7 

h 2 s 

-5.6 

-5.6 

-5.6 

HCN 

-6.4 

-6.3 

-6.4 

HNC 

- 6.6 

- 6.6 

- 6.6 

OCS 

-5.5 

-5.5 

-5.5 

HCO+ 

- 10.8 

- 11.2 

-10.9 

HCO 

- 20.8 

- 20.8 

-20.3 

HNO 

-7.9 

- 8.1 

-7.5 

CS 

- 8.0 

- 8.1 

-7.5 

HCS 

-9.6 

-10.3 

-9.2 

HCS+ 

- 12.6 

- 12.8 

- 12.1 

OCN 

- 10.2 

-10.3 

-9.9 

o 2 

- 10.8 

- 11.2 

- 11.2 

CN 

-7.9 

-7.7 

-7.2 


nuclei, and is taken to be 77 = 3.3 X 10 -12 , while N g (i) is 
the average number of molecules of the zth species adsorbed 
on to a grain. However, gas-phase chemical reactions do still 
occur for these species. For example, the dominant formation 
reaction at low temperatures for SO is O combining with 
HS. But at low temperatures (^ 30 K) SO is almost entirely 
frozen-out, which means the grain abundance dominates, 
and given the relatively large initial ice abundance (see Table 
1 ), the effect of this chemical reaction is indistinguishable in 
the total SO abundance. Likewise, CO is only frozen-out 
on to grains below 20 K and hence it is so abundant in the 
gas-phase that the effect of the shocks appears negligible. 
Indeed for species with total fractional abundances above 
approximately 10 -7 , which is most of the species injected 
into the disc at initialisation, gas-phase reactions enhanced 
by shocks are too slow to significantly affect the evolution 
of the abundances. However, although the total fractional 
abundance of some species does not depend on shocks, the 
distribution of gas-phase material is still being affected by 
spiral arms. 


3.1 Comparison with the more massive disc 

Most of the species reach the same maximum gas-phase 
abundance for both fluid elements during their respective 
shocks. This is because, despite the contribution from gas- 
phase chemical reactions, desorption dominates the forma¬ 
tion rates near the desorption temperatures (see Appendix 
A), which is exceeded for most species during the evolution 
of both fluid elements. Although the selected fluid elements 
trace the outer regions of both discs, this result is consistent 
across the entirety of both discs. Table 2 shows the max¬ 
imum gas-phase fractional abundances for the 20 species 
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Figure 10. Gas-phase fractional abundances of fluid elements extracted from the 0.17 Mq (left) and 0.39 Mq (right) discs that follow 
similar trajectories within 2.7ORPs at 30au in each disc (see Figure 4). For each disc, the upper two panels feature species whose total 
abundance, comprised of the number of molecules on grains and in the gas-phase, is constant over the duration of the simulation. 


selected at t — t 0rp , and most species show unremarkable 
differences between the lower mass and higher mass discs. 
Therefore, these results suggest that the difference in disc 
mass and instability strength does not significantly affect 
the chemistry. Rather, the chemistry is affected merely by 
the presence of shocks, and perhaps only on a short time- 
scale. 

In order to understand the features expected to exist 
in observed gravitationally unstable discs, we generated col¬ 
umn density maps for the species discussed at t — t 0rp . A 
particular aim of this section is to characterise the affect of 
system mass on chemical diversity, so these maps are com¬ 
parable to Figures 7 and 8 featured in 12011. 

The column density is defined as 

N(i,x,y) = J n(i,x,y,z)X(i,x,y,z)dz. (1) 

We obtain the number density at each x, y, z using the full 
hydrodynamics simulation as this affords us the highest res¬ 
olution possible. We then interpolate the fluid element abun¬ 
dances to the same grid using the GRIDDATA and SIMPS mod¬ 
ules in Python. To obtain results for the full vertical extent 
of such a disc, we assume that number density and abun¬ 


dances are even functions of disc height, i.e. mirrored about 
the z = 0 plane. 

Figures 11 and 12 show the column density maps for 
the selected species at t — t 0rp , separated by the signifi¬ 
cance of gas-phase chemical reactions as discussed. Figure 
11 features the species only affected by adsorption and des¬ 
orption processes due to their high initial abundances, which 
is reflected in the column density scales. We have included 
HCN, HNC and SO 2 in this figure because the most domi¬ 
nant reactions, electronic recombination of HCNH + and SO 
combining with O or OH, are only prevalent at the lowest 
temperatures and hence only have a marginal effect on the 
total abundance. Figure 12 shows the species whose abun¬ 
dances are significantly affected by chemical reactions during 
the disc evolution. 

The CO column density map is the most extensive, 
which is expected because the freeze-out temperature is 
around 20 K, and this only occurs in the outermost regions 
of the disc (> 45au). The relatively abrupt edges seen in 
this map are primarily due to the limited extent of the fluid 
elements, as the size of the full simulation is significantly 
larger than the volume defined by the fluid elements fol¬ 
lowed throughout the disc. We chose to avoid extrapolating 
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Figure 11. Logarithmic gas-phase column densities of species at t — t 0rp with gas-phase fractional abundances determined primarily 
by thermal adsorption and desorption processes. 
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Figure 12. Logarithmic gas-phase column densities of species at t = t 0rp with gas-phase fractional abundances significantly affected by 
chemical reactions. 


the abundances of species beyond the maximum radial ex¬ 
tent of the fluid parcels, as the reliability of the chemical 
results in this regime can not be assured. 

The H 2 O map is the most confined and this is because 
water possesses the highest binding energy to grains of all 
the species we consider. Hence, water requires the most en¬ 
ergy to be efficiently desorbed into the gas-phase, requiring 
a temperature greater than 150 K at these pressures. As can 
be seen in Figure 2, this only occurs for the torus of material 
approximately 10 au from the centre. 

The remaining species with the highest abundances and 
therefore highest column densities have maps that trace re¬ 
gions of the disc in between these two extremes, depen¬ 
dent on their desorption temperatures. As the highest abun¬ 
dances for most of these species exist in the hottest re¬ 
gions nearest the protostar, the edges of the column density 


maps effectively define the midplane snow lines in the xy 
plane. Hence, simulating the molecular emission from the 
disc chemistry and comparing to observations could provide 
an indicator to the reliability of our model. 

The column density maps in Figure 12, for the species 
that are significantly affected by gas-phase reactions, show 
some different features compared with the maps shown in 
Figure 11. Firstly, all of the maximum column densities are 
much lower. Secondly, CS and OCN appear to have simi¬ 
lar extent and structure to species such as CO 2 and SO, 
however, their maximum column density occurs in a torus 
at approximately 20 au from the centre. Furthermore, the 
HCO + map spans the entire fluid element distribution but 
possesses a large inner hole. This is due to HCO + being pri¬ 
marily formed through a reaction with CO, but destroyed 
through charge exchange with H 2 O and reactions with HCN 
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and HNC, all of which are most abundant in the inner disc 
region. HCS + is also removed by charge exchange with H 2 0 
but reflects the distribution of CS and hence is not as ex¬ 
tensive. Aside from HCO + , CN is one of the only species 
to show clear evidence of spiral structure. While not a com¬ 
ponent of our initial chemical abundances, CN is initially 
formed in the gas phase via cosmic-ray induced photodis¬ 
sociation of HCN and HNC (see, e.g., Gredel et al. 1989). 
The desorption temperature of CN is around 120 K, which 
occurs in these spiral features due to the shock heating. At 
higher temperatures the reaction CN + NH 3 —HCN + NH 2 
destroys the molecule, hence producing a very large inner 
hole. As CN maps areas of moderate temperatures, it could 
be an ideal tracer for transient heating events in protoplan¬ 
etary discs, potentially driven by gravitational instabilities. 

Although the less massive and more massive discs rep¬ 
resent different strengths of GI activity, by using a compa¬ 
rable dynamical age in both systems, defined as t — t orp 
(2.7ORPs at 30au), we are able to draw some comparisons. 
Figures 7 and 8 in 12011 show column density maps for the 
same species we have included in this work, except for CN 
and C0 2 . All of the maps included in the 12011 paper show 
much more defined spiral structure, but this is due to the 
more massive system driving stronger density waves than in 
the less massive disc. Despite this there are very few remark¬ 
able differences between the sets of column density maps; 
the proportional extent of species in relation to the total 
disc size is comparable and the maximum column densities 
are largely coincident. For example, the CO column density 
maps show much less spiral structure in the lower mass disc 
at t — t 0 rp, but this is purely a consequence of the mass and 
instability strength differences between the discs; the max¬ 
imum CO abundance is identical in both discs due to the 
consistent initial conditions and low freeze-out temperature 
(see Table 2). There are, however, some distinct differences 
in the column density maps of HCO + . 

There is a 2.5 times lower peak HCO + abundance in the 
lower mass disc, which also manifests itself in the column 
density maps. This is a consequence of the stronger spiral 
density waves and larger radius of the more massive disc. 
HCO+ is formed through a reaction involving CO, which 
proceeds at the fastest rate in low to moderate tempera¬ 
tures (below approximately 150 K). The destructive reac¬ 
tions HCO+ + HNC ->• HCN + H + + CO, and more notably, 
HCO+ + HCN -► HCN + H + + CO become significant at 
temperatures around 40-60 K. In the outer regions of the 
lower mass disc these reactions are prevalent due to the 
weak shock strength. However, the shocks are stronger in 
the higher mass disc, and hence, in the outer regions, these 
destructive reactions are much less significant. This results 
in a higher peak HCO + abundance. 

Moreover, the different masses of the discs appear to 
affect the size of the HCO + inner hole. At higher temper¬ 
atures, occurring in the inner disc, H 2 0 destroys HCO + at 
a dominant rate, producing an inner hole depleted of the 
cation. After 2.7 ORPs, the maximum temperature in the 
more massive disc is 15 per cent hotter than in the lower 
mass disc and the extent of the high-temperature region in 
the inner disc is larger. Therefore, in the lower mass disc, 
HCO + is suppressed significantly at radii approximately 1.4 
times smaller than in the massive disc (see Appendix B, 


Figure Bl). This result suggests that HCO + could poten¬ 
tially be used to characterise disc mass and the strength of 
instabilities in young systems. 

The most striking difference, in terms of tracing spi¬ 
ral structure, between our column density maps and those 
featured in 12011 occurs for HCS. HCS only traces the in¬ 
nermost disc regions in the lower mass disc, whereas, in 
the more massive disc, HCS also maps the strongest spiral 
features. More importantly, however, the HCS abundance 
in these strongly shocked, outer disc regions is equivalent 
to the peak abundance in the innermost disc; a significant 
number of other species have distinctly lower abundances in 
the outer spiral structure. This is because HCS is desorbed 
above 65 K (see Table Al), but once in the gas-phase, is not 
destroyed at significant rates. Hence, gas-phase HCS is pre¬ 
served in shocked regions. As a result, HCS could be a useful 
tracer of spiral shocks in particularly massive systems. 


3.2 Chemical evolution of the lower mass disc 

By simulating the 0.17 M® protoplanetary disc for multi¬ 
ple dynamical time-scales we are able to investigate how 
a young, gravitationally unstable system, potentially analo¬ 
gous to the early Solar System, could develop. In addition to 
the changes in physical properties, such as enhanced maxi¬ 
mum density and temperatures as discussed, it is interesting 
to observe the effects on the chemical evolution of the disc, 
and in particular characterise whether these are short-term 
results or can leave a lasting imprint on further chemical 
evolution. 

Figure 13 shows the evolution of the fractional gas- 
phase abundances of the same 20 species used in Section 
3.1 for fluid elements in the outer and inner regions of the 
lower mass disc. The fluid elements are the same as featured 
in Figures 8 and 9 and orbit at approximately 30 au and 
6 au, respectively. The effects of the shocks on the outer disc 
fluid element are clearly seen throughout the entire dura¬ 
tion of the abundance histories, except for CO due to its 
very low freeze-out temperature. The shock heating raises 
the temperature of the fluid element, affording more energy 
for species to desorb or overcome activation barriers. The 
exception to this trend is HCO + , which as discussed is de¬ 
stroyed by high abundances of H 2 0, HCN and HNC. 

The abundance histories of the inner disc fluid element 
appear significantly different but are entirely consistent with 
the outer disc abundances. Every species effectively reaches 
its maximum abundance once the temperature exceeds the 
energy required for desorption. This occurs at different tem¬ 
peratures, and hence at different times, for each species. 
For example, HCN and HNC are efficiently desorbed above 
60 K, which is instantly reached upon initialisation for the 
inner disc fluid element. Hence, the fractional gas-phase 
abundance of HCN and HNC is instantaneously saturated 
and persists until the end of the simulation. H 2 0, on the 
other hand, is desorbed above 150 K, which is only consis¬ 
tently exceeded after around 600 yr. Before this time, the 
shocks can be seen to have an effect on the abundance as 
already discussed, and afterwards effectively all the water 
is in the gas-phase. The behaviour of CN is more complex, 
but has already been explained. CN is efficiently desorbed 
above 45 K and hence reaches maximum gas-phase abun- 
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Figure 13. Gas-phase fractional abundances of fluid elements extracted from the inner (left) and outer (right) regions of the 0.17M© 
disc (see Figure 8). For each region, the upper two panels feature species whose total abundance, comprised of the number of molecules 
on grains and in the gas-phase, is constant over the duration of the simulation. 


dance nearly instantly. However, at high temperatures CN 
is formed through cosmic-ray induced photodissociation and 
destroyed by NH3. Hence, once virtually all of the CN is in 
the gas-phase, chemical reactions occur at a rate significant 
enough to affect the abundance. There are more apprecia¬ 
ble differences in species for which we have not displayed 
fractional abundances, but we choose to exclude these as 
they have not been observed in protoplanetary discs and/or 
have an insignificantly low abundance. For example, OH has 
been detected in protoplanetary discs (e.g. Fedele et al. 2012; 
Meeus et al. 2012; Fedele et al. 2013) but we only recover a 
peak OH abundance of 10“ 14 because we assume our disc is 
well shielded from UV photons, which corresponds to a peak 
column density of 10 12 cm -2 that is likely undetectable. Fur¬ 
thermore, although reactions of atomic hydrogen can play 
a destructive role at temperatures greater than 300 K, the 
fractional abundance of atomic hydrogen in the gas phase 
is approximately constant throughout the disc evolution at 
~ 10 —11 , which corresponds to approximately lppcm -3 . As 
such, this relatively small reservoir of H is quickly used up 
and is not available for further reactions. Other reactions 
involving O and C atoms play a more dominant role in the 
chemical evolution. 

In order to further assess the implications of GIs on the 


chemical evolution of young, self-gravitating protoplanetary 
discs, we produced column density maps at t — for com¬ 
parison with Figures 11 and 12. We used equation 1 and the 
same procedure and presentation as discussed in Section 3.1. 
Every species shown in Figures 14 and 15 has a larger max¬ 
imum column density compared to t — t 0rp , which can be 
attributed to a combination of enhanced maximum nuclei 
density and maximum temperature. This suggests, there¬ 
fore, that despite the fact that some species are destroyed 
in the inner disc, we still expect to detect higher column 
densities of all species in a hotter and denser disc. 

Investigating the structure of the column density maps 
at t — £fi n we see that CO is again the most extensive species 
because of its highly volatile nature. However, due to the 
more flocculent structure evident in the nuclei column den¬ 
sity map at the end of the simulation (see Figure 6), the 
spiral structure is much more prominent. This strengthening 
of spiral features is also observed in species with relatively 
low desorption temperatures, such as CO 2 , H 2 CO and CS, 
but species with the highest binding energy to grains, such 
as H 2 O and NH3 still only trace the innermost regions of 
the disc. CS in particular appears to be a useful molecule in 
tracing the spiral features expected in young, self-gravitating 
discs but only after the spiral waves have fully developed. 


MNRAS 000 , 000-000 (0000) 


































































































































GIs in a protosolar-like disc I 13 


log N [cm- 2 ] 




0 


20 


: 19 


-60 -40 -20 0 20 40 60 

x [au] 



D 


21 


20 


D 


23 




: 20 


Figure 14. Logarithmic gas-phase column densities of species at t = £fi n with gas-phase fractional abundances determined primarily by 
thermal adsorption and desorption processes. 
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Figure 15. Logarithmic gas-phase column densities of species at t = £fi n with gas-phase fractional abundances significantly affected by 
chemical reactions. 


CN, on the other hand, could prove to be excellent in char¬ 
acterising instabilities in embedded discs as it traces non- 
axisymmetric features throughout the entire disc evolution 
with a relatively high column density, and, more importantly 
for current observational constraints, in the outer regions of 
the disc. 

3.2.1 Persistent effects of shock heating 

The shock heating in the disc clearly has an instantaneous 
effect on the chemistry, but it is worthwhile to investigate 
whether these effects persist. In the inner disc, once the des¬ 
orption temperature is exceeded for a particular species, the 
abundance appears unaffected by further shocks. However, 
Table 2 shows that some species have enhanced maximum 
abundances at the end of the simulation, e.g. about a 10 


times increase in maximum abundance of HCS, despite the 
high temperatures and densities of the inner disc at t — £ 0 r P - 
This suggests that the enhanced densities and temperatures 
of the inner disc regions driven by GIs in a young, massive 
disc can lead to subtle changes in the chemical composition. 
The inner regions of discs (i.e. within 20 au), however, par¬ 
ticularly in young, embedded systems, are very challenging 
to observe. 

In the more observationally accessible outer disc, the 
shocks have a continual effect on the chemical abundances 
due to the low ambient temperature. Although most species 
appear to reach the same abundances post-shock as pre¬ 
shock, we investigated further and found evidence to suggest 
that consistent shock-heating changes the chemical compo¬ 
sition, at least over the short timescale of our simulation. To 
do this we took a fluid element in the outer disc and fitted 
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Table 3. Fractional gas-phase abundances at t = 2000 yr for a 
fluid element in the outer disc with shocks (WS), and for the same 
fluid element without shocks via fitting a low-order (WOSLO) 
and high-order (WOSHO) polynomial to the temperature and 
number density minima. Note that we present the results for t = 
2000 yr because at the end of the simulation the fluid element is 
encountering a shock. 


Species 

i 

WS 

log X (i) 
WOSLO 

WOSHO 

CO 

-4.4 

-4.4 

-4.4 

C02 

-19.1 

-17.3 

-17.4 

SO 

-21.7 

-19.7 

-19.7 

S02 

-27.3 

-24.5 

-24.4 

NH3 

-20.3 

-24.8 

-24.8 

H20 

-17.5 

-16.2 

-16.2 

H2CO 

-12.9 

-11.9 

-11.9 

H2S 

-23.1 

-21.3 

-21.4 

HCN 

-23.8 

-24.1 

-24.1 

HNC 

-23.7 

-24.8 

-24.9 

OCS 

-21.3 

-20.1 

-20.1 

HCO+ 

-11.8 

-11.5 

-11.5 

HCO 

-26.5 

-25.0 

-25.1 

HNO 

-21.2 

-24.3 

-24.2 

CS 

-22.5 

-20.9 

-21.0 

HCS 

-22.1 

-20.8 

-20.9 

HCS+ 

-23.6 

-22.2 

-22.3 

OCN 

-27.8 

-28.3 

-28.3 

02 

-18.0 

-15.2 

-15.2 

CN 

-18.0 

-21.3 

-21.3 


a low and high-order polynomial to binned minima of the 
temperature and number density (see Figure Cl). This al¬ 
lowed us to investigate the chemical evolution of an artificial 
fluid element more appropriate for a trajectory within a disc 
absent of shocks in the outer regions. 

We found that the abundances of most species are not 
affected significantly when shocks are removed, as the abun¬ 
dances near the end of the simulation are comparable, as 
shown in Table 3. Hence, the chemical changes for most 
species are not persistent. However, HNO, CN and NH3 
have significantly higher abundances near the end of the 
simulation when the fluid element is consistently shock- 
heated in the outer disc. For example, devoid of shocks, 
HNO is formed via O + NH 2 —t HNO + H at an approx¬ 
imate rate 10 -31 cm -3 s -1 . In shock-heated regions, how¬ 
ever, this rate reaches a maximum of 10 -22 cm _3 s _1 , and 
the frequency of shocks results in a limited minimum rate of 
lO -28 cm _3 s _1 . This leads to a HNO abundance enhanced 
by a factor of 10 3 . NH3 is enhanced by a factor of 3 x 10 4 
due to this same phenomenon, only the rate of the reaction 
NH 4 + + e _ —»• NH3 + H is consistently enhanced at least 
10 4 times from shock-heating. CN is enhanced by a factor of 
2 x 10 3 in much the same manner, except that desorption is 
the dominant formation reaction due to the lower desorption 
temperature of CN (see Table Al). 

Shocks in a relatively low-mass, protosolar-like disc are 
seen to cause persistent changes in the abundance of some 
species. This has several implications for future studies of 
young protoplanet ary discs. Firstly, whether this effect is 
independent of the initial conditions should be investigated, 
because subtle differences in early-phase composition may 


not be preserved through a short evolutionary period, re¬ 
sulting in similar abundances at the end of the GI phase 
across a variety of young systems. Secondly, if gravitational 
instabilities are a significant phase of early disc evolution, 
simulations of more evolved systems should consider the 
processing of GIs before setting their initial abundances, 
in order to account for the enhancement in some impor¬ 
tant species. Ideally, models of protoplanetary disc chem¬ 
istry should use a self-consistent approach concerning the 
molecular core collapse phase (e.g. Visser et al. 2011), fol¬ 
lowed by Gl-driven evolution if the disc-to-star mass ratio 
is large enough. The resulting abundances can then be used 
as initial conditions for the simulation of a more evolved, 
low-viscosity or passive disc, at which point an axisymmet- 
ric <a-disc prescription becomes appropriate. Furthermore, 
these results suggest that observations of discs could still 
be used to characterise disc dynamics even if they can not 
be resolved spatially, which is an important finding. This is 
because successive shocks enhance the abundances of some 
important species, and so an observation of a species such as 
CN, that may otherwise be undetectable in the outer regions 
of a more passive disc, indicates an enhancement in temper¬ 
ature that could be produced by GIs. Moreover, even if the 
effects of shock heating are short-lived for some species, the 
enhancement in abundance is very significant, e.g. up to a 
factor of 10 20 for SO 2 (see Figure 13), Therefore, detections 
of transient species may still be indicative of GIs. Further 
investigation is warranted in this area to determine how per¬ 
manently gravitational instabilities affect disc composition 
over longer timescales. 

3.3 Comparison with axisymmetric models 

Comparing the results reported in this paper to other results 
in the literature is not straightforward due to the varying 
simulation dynamics, techniques and assumptions. We have 
used the chemical network developed in 12011, yet even a 
comparison with this work is complicated because the less 
massive and more massive discs represent different stages 
of GI activity and occupy different parameter spaces. Nev¬ 
ertheless, Table 4 shows maximum abundances taken from 
two recent studies, Walsh et al. (2010) (hereafter W2010) 
and Akimkin et al. (2013) (hereafter A2013), to provide a 
very rough comparison to the maximum abundances at the 
end of our lower mass disc simulation. 

We extracted the peak abundances from W2010 by vi¬ 
sually inspecting their ‘fiducial’ model between 10 au and 
50 au, which introduces some uncertainty into the quoted 
values. Nevertheless, there are significant differences be¬ 
tween our results that arise from a multitude of factors. 
Firstly, W2010 model a 1 Myr old laminar model disc, which 
is more evolved than the system we consider and features 
only local mass transport. As a result their disc possesses the 
vertically inverted temperature structure observed in Class 
II YSOs due to the photon-dominated region near the disc 
surface. Therefore the highest abundances for most species 
in W2010 occur in a ‘transition’ layer confined to z/r « 
0.3, where UV and X-ray dissociation is pertinent, which 
explains their greatly enhanced OH abundance. If we com¬ 
pare only within the midplane regions we recover a simi¬ 
lar peak OH abundance. This is also true for the HCO + 
abundance, which is expected to be due to the ubiquity 
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Table 4. Maximum gas-phase abundances reported within our 
work, W2010 and A2013. A dash represents a missing reported 
value for that species in the corresponding work. 


Species Maximum log X(i) 


i 

This work 

W2010 

A2013 

CO 

-4.4 

-4 

-4.0 

C02 

-4.4 

-6 

- 

NH3 

-5.5 

- 

-4.7 

H20 

-3.7 

-4 

-5.0 

H2CO 

-5.7 

-9 

-9.0 

HCN 

-6.4 

-7 

-8.5 

HCO+ 

-10.9 

-6 

- 

cs 

-7.5 

-8 

- 

02 

-11.2 

- 

-4.5 

CN 

-7.2 

-7 

- 

C2H 

-9.9 

-7 

- 

N2H+ 

-18.9 

-11 

-11.0 

OH 

-13.9 

-4 

- 


of gas-phase CO in both models. The midplane abundance 
of H 2 CO in W2010 is much lower than we report due to 
the more efficient freeze-out in their colder midplane and 
lower initial abundance; W2010 use significantly different 
initial abundances, focusing on oxygen-rich low-metallicity 
elements rather than comet ary ices. Due to the strong de¬ 
pendence on initial abundances for some species, those that 
are only significantly affected by adsorption and desorption 
processes, we find a higher abundance at the end of our sim¬ 
ulation. Furthermore, we do not include atomic N initially in 
our chemical network which results in a significantly lower 
N 2 H+ than W2010 found. 

The maximum abundance values quoted for A2013 were 
determined from figures showing the abundances at 10 au 
and 50 au within their disc model. We use different initial 
abundances in our model than A2013, who use oxygen-rich 
low-metallicity abundances similar to W2010, which helps 
explain our enhanced values for H 2 O, H 2 CO and HCN. Con¬ 
trary to A2013, we do not include N or N 2 in our initial con¬ 
ditions and furthermore, A2013 incorporate grain growth in 
their simulation, which they report enhances the abundance 
of nitrogen atoms in the midplane. Hence, we record a lower 
N 2 H + abundance. Interestingly, however, A2013 report the 
same peak N 2 H + abundance as W2010 despite the lack of 
grain growth treatment in the latter. Our lower maximum 
abundance of O 2 is a consequence of our model contain¬ 
ing 1.3 times less elemental oxygen and the A2013 model 
including UV and X-rays that enhance the abundance of 
OH in the upper disc, which subsequently leads to synthe¬ 
sis of molecular oxygen via a neutral-neutral reaction with 
atomic oxygen. As the A2013 disc follows the classic a pre¬ 
scription with an array of photoprocesses incorporated into 
their chemical model, they recover a heated layer in the disc 
possessing the richest chemistry, similar to W2010, which is 
common of more passive discs. 

Realistically, these comparisons are very weak be¬ 
cause of the significantly different dynamics between our 
model and other models in the literature. In essence, what 
they show is that in a more passive, axisymmetric <a-disc, 
a photon-dominated, chemically rich layer sits above a 
shielded, cold, chemically sparse midplane. However, in a 


non-axisymmetric, gravitationally unstable disc, global in¬ 
stabilities and mass transport produce a hot midplane that 
experiences a rich chemical evolution. Ideally, other 3D sim¬ 
ulations of Gl-driven discs need to be completed before ap¬ 
propriate comparisons can be made to the results presented 
in this paper. 


4 CONCLUSION 

We have modelled the physical and chemical evolution of a 
0.17 Mq protoplanet ary disc over a period of approximately 
2000 yr. The disc surrounds a 0.8 M 0 protostar that will 
evolve into a Solar-like star. As a result, our work may be 
indicative of an early dynamical and chemical phase of evo¬ 
lution of our Solar System. 

This paper extends previous work that suggested grav¬ 
itational instabilities significantly affect the chemistry in 
more massive systems. The main results that we have found 
from simulating a lower mass, protosolar disc are as follows: 

• Spiral waves generated in the disc heat material 
through shocks that enhance the local temperature and in¬ 
crease the rates of desorption and some endothermic reac¬ 
tions. This subsequently increases the gas-phase abundance 
of most species across the entire disc, including the mid¬ 
plane. 

• In the more massive disc, the majority of species 
trace distinct spiral structure in the column density maps 
after 2.7 ORPs at 30 au (defined as t — t orp )• After the same 
number of outer rotational periods in the lower mass disc, 
the majority of species do not trace much spiral structure 
because the spiral density waves are of lower amplitude. 

• At the end of our simulation spiral features driven 
by GIs in a protosolar disc are traced by a significant num¬ 
ber of important species; CN is a particularly strong tracer 
as it is most abundant in the outer disc regions and has a 
relatively high column density. 

• Comparing the less massive and more massive discs 
at t — torp reveals that most species have reduced or en¬ 
hanced peak abundances by a factor of only two to four 
times. Therefore, this suggests that the presence of gravi¬ 
tational instabilities in young, massive protoplanetary discs 
defines the chemical evolution more significantly than the 
mass differences of the discs. 

• We have found that GIs, due to global transport 
of disc material, produce a chemistry-rich midplane. This is 
in contrast to studies of more evolved, axisymmetric discs, 
that report peak abundances in a heated layer above the 
cold ‘gas-phase desert’ midplane, caused by photoprocesses 
and local transport of material. 

• Rapid succession of shocks can lead to long-lasting 
changes in the inner disc composition. At small radii, the 
shock heating quickly increases the temperature beyond the 
desorption energy for all neutral species. Therefore we re¬ 
cover abundances close to the protostar that are likely higher 
than in a more passive disc. 

• In the outer disc, the effect of shock heating is easily 
distinguishable because of the lower ambient temperature. 
For most species this effect is short-lived, however, successive 
shock heating permanently alters the abundances of HNO, 
CN and NH 3 in the outer disc over the duration of our sim- 
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ulation. This has two major implications. Firstly, chemical 
simulations of more evolved systems should perhaps con¬ 
sider the processing of GIs when setting initial conditions, 
if indeed GIs are present in the early phases of disc evolu¬ 
tion. Secondly, observations of discs could still be used to 
characterise disc dynamics even if they can not be resolved 
spatially as shocks cause transient and permanent changes 
in gas-phase abundances. 

• H 2 O destroys HCO + via charge exchange interac¬ 
tion, creating a hole in HCO + in the inner disc. The extent 
of H 2 O is predominantly determined by instability strength, 
which is characterised by disc mass. Therefore, observations 
and measurements of the HCO + inner hole in real systems 
could be used as an additional method for characterising 
protoplanet ary disc mass. 

We have not incorporated envelope accretion or out¬ 
flows into our model in order to focus on the effects of gravi¬ 
tational instabilities on the chemical evolution of protoplan¬ 
etary discs. This is appropriate as collimated outflows are 
assumed to have an insignificant effect on disc composition 
and a surrounding envelope would likely only enhance the 
molecular abundances in the outermost regions of the disc. 
Sakai et al. (2014) have found an enhancement of SO at the 
centrifugal barrier of a Class 0 object that could be produced 
from a weakly shocked region as envelope material accretes 
on to the disc whilst exceeding the Keplerian velocity. Our 
results support this theory as the gas-phase abundances of 
volatile species are enhanced through shock heating. How¬ 
ever, as the authors find no significant SO abundance in¬ 
terior to the centrifugal barrier, the impact from this phe¬ 
nomenon on our results concerning gravitational instabilities 
is likely negligible. 

Recent observational advancements such as ALMA will 
revolutionise the studies of protoplanetary discs, and open 
up the possibilities for investigating young, embedded ob¬ 
jects. The resolutions attainable suggest that if the 0.17 M 0 
disc was located within the Taurus-Auriga cloud complex, 
its spiral structure should be distinguishable, which is an ex¬ 
citing prospect for characterising gravitational instabilities 
in real systems. As a cautionary note, however, it is possible 
to have well defined spiral structure and not detect it, or for 
complicated spiral structure to give an incorrect appearance 
at low resolution. For example, CO appears to be the best 
species for directly observing spiral structure in real systems, 
but at low resolution the image could appear as two rings, 
which would be incorrectly interpreted as a gap. Therefore, 
in paper II of this series, we will perform radiative transfer 
calculations in order to assess the detectability and image 
fidelity of spiral structure in a lower mass disc, focusing on a 
combination of the species we have found to be encouraging 
tracers such as CN, HCS and CO. We will then compare 
these to the results of Douglas et al. (2013) in an attempt to 
determine the importance of system mass on observability 
of such features. 
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Table Al. Desorption temperatures. 


Species 

Binding Energy [K] 

T des [K] 

CO 

855 

26 

C02 

2575 

74 

SO 

2600 

72 

S02 

3405 

93 

NH3 

5534 

151 

H20 

5773 

167 

H2CO 

2050 

57 

H2S 

2743 

76 

HCN 

2050 

57 

HNC 

2050 

57 

OCS 

2888 

79 

HCO 

1600 

44 

HNO 

2050 

57 

CS 

1900 

53 

HCS 

2350 

65 

OCN 

2400 

66 

02 

1200 

33 

CN 

1600 

44 


APPENDIX A: DESORPTION 
TEMPERATURES 


The rate of desorption from a surface is given by the Polanyi- 
Wigner relation (Tielens 2005), 

Ri = Vie~ Ei,kTsi (Al) 


where Ei is t he adsorption bind ing energy of species i, Vi 
= 1.6xl0 n ^/(Ei/k)/(mi/mu) s _1 is the vibrational fre¬ 
quency of the species in the surface potential well, k is Boltz¬ 
mann’s constant, and rrn and run are mass of species i and 
hydrogen, respectively. T gr is the temperature of the grains, 
which we assume to be in equilibrium with the gas tem¬ 
perature. The freeze-out temperature can be determined by 
equating the flux of desorbing molecules to the flux of gas 
molecules adsorbing on to the grain surface, calculated from 
gas kinetic theory. Thus, 


N s ,iRif s ,i = 0.25 mvi (A2) 


where N s ,i is the number of adsorption sites per cm 2 , f Sj i is 
the fraction of the surface adsorption sites that are occupied 
by species i (see Ilee2011, Equation 9), rii is the gas-phase 
number density of species i 1 and v\ is the thermal speed. A 
sticking probability of unity is assumed. By solving for the 
dust temperature we derive the desorption temperature for 
species i, (Hollenbach et al. 2009), 


rp ^ Ei 

J- de .s — 


In 


{ 4A T s ,jf s ,jVj \ 
V UiVi ) 


(A3) 


Using equation A3 we derived the desorption temperatures 
for the neutral species featured in this paper, which are 
shown in Table Al. We have used an average gas-phase 
number density for each species in order to provide desorp¬ 
tion temperatures that are applicable across the whole disc, 
within a maximum error of 10 per cent. 


higher mass discs at t — t 0rp . In the more massive disc the 
overall extent of significantly abundant H 2 O is much larger 
due to the higher amplitude spiral waves. As a result, the 
H 2 O and HCO + abundances are enhanced and suppressed, 
respectively, over a larger fraction of the disc radius, and 
so the inner hole seen when viewing from above is larger 
than in the lower mass disc. This result indicates that the 
measurement of the HCO + hole size in real systems could 
potentially be used as a tracer for disc mass. 


APPENDIX C: REMOVING SHOCKS 

Figure Cl shows the temperature and number density his¬ 
tory of the same fluid element in the outer disc with and 
without shocks. Artificial fluid elements were constructed 
by fitting low and high-order polynomials to the minima of 
temperature and number density, which are more appropri¬ 
ate for trajectories within a passive disc. The abundances 
of significant species were calculated for the fluid elements 
with and without shocks, and compared to determine the 
persistent effects of shock heating within a gravitationally 
unstable disc. 


APPENDIX B: HCO+ INNER HOLE 

Figure B1 shows the abundances of H 2 O and HCO + as 
a function of midplane radius within the lower mass and 
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M = 0.17M(7) 




Figure Bl. Radial extent of H 2 O and HCO+ fractional gas-phase abundances in the midplane of the 0.17M© (top) and 0.39 Mq 
( bottom) discs at t = t Qrp . The abundances for every midplane cell are plotted, and because of the shock heating and non-axisymmetric 
nature of the discs, there is a wide range in abundance values at each radii. 



Figure Cl. Temperature and number density history of a fluid 
element in the outer disc. Overplotted are polynomials of low 
and high-order fitted to the minima of temperature and density, 
representing the same fluid element history devoid of shocks. 
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